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1.0  BROADBAND  PROCESSING  ISSUES 


In  this  section,  we  introduce  issues  related  to  the  analysis  of  broadband  processes  in 
the  context  of  the  source  localization  problem  via  spatial  sampling.  We  will  present  a 
linear  systems  model  of  the  propagation  environment  and  show  how  it  influences  the 
spatial/spectral  estimation  of  these  processes.  We  also  show  how  this  linear  systems  view 
lends  itself  to  the  analysis  of  the  moving  source  problem. 

1.1  NOTATION 

In  this  document,  bold  lower  case  and  upper  case  letters  denote  vectors  and  matrices 
respectively,  x  represents  a  vertical  vector,  x*1  its  complex  conjugate  or  Hermitian  trans¬ 
pose,  and  xT  its  transpose. 

1.2  LINEAR  SYSTEMS  MODEL  OF  PROPAGATION  ENVIRONMENT 

Assume  that  an  acoustic  source  emits  a  real,  zero  mean,  wide  sense  stationary  (WSS) 
ergodic  random  process  s(t)  with  bandwidth  B.  The  process  propagates  through  some 
acoustic  environment,  and  we  spatially  sample  the  field  with  an  array  of  M  sensors.  Each 
sensor  is  indexed  by  its  location  vector  fisi  =  (rsi,zsi,0si)  where  rsi,zsi,  and  QSI  denote 
range,  depth,  and  azimuth  with  respect  to  the  source.  x(t;fisi)  represents  the  process 
received  by  the  ith  sensor  due  to  s(t). 

1.2.1  Time-Invariant  Channel  Filter  Model 

It  is  instructive  to  view  x(t\psi)  as  resulting  from  a  linear  filtering  operation  on  s(t) 
(Knight,  Pridham,  and  Kay,  1981).  Figure  1-1  represents  the  propagation  environment  as 
a  bank  of  linear  filters  where  h(f,psi )  is  the  time  invariant  channel  impulse  response  for  the 
ilh  source/sensor-pair  filter  and  y((t)  is  the  output  of  the  sensor  due  to  si)  and  additive 
noise  n^t),  which  is  also  WSS  ergodic  and  statistically  uncorrelated  with  x (r;/U).  Natu¬ 
rally,  these  “spatial  filters”  are  causal  and  stable.  We  also  will  assume  that  the  channel 
band  widths  are  much  greater  than  B.  Figure  1-2  depicts  a  sample  time  invariant  impulse 
response.  We  can  imagine  sending  an  impulse  down  a  channel  and  receiving  this 
response.  Here  the  time  spread  L,  -  £>.  generally  can  be  interpreted  as  the  result  of  multi- 
path  (Kennedy  and  Lebow,  1964). 

If  o]( w)  is  the  power  spectral  density  of  the  process  s(t),  we  know  that  the  power 
spectral  density  of  x(t,psi)  is 


CTW^<(w)  =  I//(®;0h-)I2<7?(®)  (1.1) 

and  the  cross  power  spectral  density  of  and  x(t;fisk)  is 
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ahsixPsSa^  =  H(<°'’Psi)H  *(<*>■,  Psk)oK<»)  (1.2) 

•th  th 

where  H(oi\f}Si)  and  H(u\f}sk)  are  the  frequency  responses  of  the  i  and  k  source/ 
sensor-pair  filters  respectively. 

1.2. 1.1  Array  Processing 

We  first  will  be  concerned  with  frequency-domain  array  processing  applications 
(Dudgeon  and  Mersereau,  1984,  and  Gingras,  1989)  since  most  well-developed  acoustic 
propagation  codes  directly  compute  the  channel  frequency  responses  rather  than  impulse 
responses.  We  estimate  the  temporal/spatial  spectrum  />(o);/?)  of  the  field  by  computing 
the  quadratic 

^con  (1-3) 

for  the  conventional  case  and 


^mvdr  P)  ~ 


for  hie  minimum  variance  distortionless  response  (MVDR)  processor.  Here  e^-p)  is  the 
steerir  g  vector  where  the  ith  element  is  an  estimate  of  the  channel  frequency  response 
sampled  at  w  for  the  ith  sensor  indexed  by  some  spatial  vector  f3,  (indicating  the  “look” 
direction)  and  /?(w)  is  the  cross  spectral  density  matrix  (CSDM),  defined  by 


*(<u)  =  lim  \ ■E{y(o})yH{a)))  (1.5) 

T  — ►  oo  l 

where  the  i'h  element  of  y(w)  is  given  by  the  stochastic  integral 

77  2 

y,(v)  =  J  yi(t)e'jMdt  (1.6) 

-7/2 

and  E{  }  denotes  the  expectation  operator.  For  the  case  of  figure  1-1  and  using  equations 
1.1  and  1.2,  /?(w)  can  be  written  as 


R((o)  =  o25(a))h{0),p5)hH{0}.ps)  +  Q( (o)  .  (1.7) 

We  refer  to  h(o>;ps)  as  the  signal  vector  where  the  i,h  element  is  and  Q(o>) 

represents  the  CSDM  of  the  noise  process  n;(t).  We  compute  equations  1.3  or  1.4  for 

various  (w;/?),  and  if  e^-p)  = 

1.2. 1.2  Plane-Wave  Case 

The  above  analysis  also  applies  to  the  well-known  plane-wave-based  array  processing 
problem.  In  this  case,  the  impulse  responses  of  figure  1-1  will  be  delta  functions  shifted  in 


^ we  will  get  a  peak  on  the  surface. 

h(0);ps) 
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time  t  by  an  amount  corresponding  to  the  direct  path  delay.  For  the  case  in  figure  1-2, 
h(t;0i )  would  be  equal  to  8 (t  -  D)  and  the  elements  of  h((jr,Ps)  would  be  the  corresponding 
phase  shifts. 

1.2. 1.3  Estimating  the  CSDM 

We  know  that  since  s(t )  is  WSS  and  ergodic  and  is  stable  and  time  invariant, 

■*(t  \fisi)  and  y(f()  also  are  WSS  ergodic.  So  the  ensemble  average  of  equation  1.5  can  be 
estimated  using  frequency  averages  by  dividing  the  sensor  data  into  K  sufficiently  long 
segments  and  computing 

R(«>)  =  (1-8) 

K  i 

th 

where  y^w)  is  calculated  from  the  normalized  Fourier  transform  of  the  /  segment  of 
data  (Johnson,  1982).  Generally,  as  K  and  the  respective  transform  length  become  larger 
and  the  time  gap  between  the  Ith  and  the  (/  -  l)(/!  Fourier  transforms  becomes  smaller, 
R( to)  becomes  a  “better”  estimate  of  tf((o).  Usually  K  >  M. 

1.2. 1.4  Fourier  Transform  Length/Spectral  Sampling  Rate 

The  method  of  estimating  the  CSDM  from  equation  1.8  requires  us  to  compute 
Fourier  transforms  that  are  “sufficiently  long.”  Aside  from  the  statistical  requirements 
needed  to  estimate  the  source  power  spectral  density  (oj(to)  in  equation  1.7),  the  neces¬ 
sary  length  must  be  a  function  of  the  array  geometry  and  the  propagation  environment. 
To  see  this  let  s(t)  be  a  broadband  process  and  let 


''  +  r/  ^ 


xti((o)  =  -Jff  \  \s(t-T)h(r,fis)e~J0>ldtdt 


'  J  t ,  Di 


Li  It+Tf 

-j=  j  h(r ;  psi)  j  s(t  -  r)  e  J0>,dr  dt 

*  f  Dt  t, 


(1.9) 


be  the  signal  component  of  the  ilh  element  of  the  /"’vector  from  equation  1.8.  Now,  if 


S,(co)  = 


(1.10) 


then 

xu((o)  =  Si{(o)  ^h{x\P5)  e'^dx  (i.n) 

D, 

if  Tj  »  Lr  Then 
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(1.12) 


xu(o>)  *  Sfa)  H((o;Ps)  . 

Now,  since  we  really  are  interested  in  correlations  between  sensors,  i.e. ,  phase  differ¬ 
ences  between  elements  of  x /  (a>)  as  in  equation  1.2,  the  condition  that  Tf  »  Li  might  be 
too  strict.  We  can  imagine  time-shifting  the  impulse  responses  of  the  source/sensor  pair 
filters  by  an  amount  Ds=  min (Z),,D2,...,DM  ).  Then  if  Ls=  ma x(LvL2,...,LM),  we  can 
approximate  good  phase  differences  for  all  sensor  pair  elements  of  the  CSDM  if 


T»L  -D  . 

f  s  s 


(1.13) 


This  result  is  analogous  to  the  plane-wave  case  where  the  length  of  the  Fourier  transform 
is  required  to  be  much  longer  than  the  maximum  propagation  time  across  the  array  (Wax, 
Shan,  and  Kailath,  1984). 

Interestingly,  this  requirement  is  less  important  for  narrowband  processes  (Wax,  Shan, 
and  Kailath,  1984)  and,  in  fact,  is  unnecessary  for  sinusoidal  processes.  To  see  this, 
assume  s(t)  is  a  complex  sinusoidal  process  where  0  is  some  random  phase.  Then 

equation  1.10  becomes 


(1.14) 


where  c  is  a  complex  constant.  Then  equation  1.11  becomes 

xu(w)  =  S,((0)  j  h(r,fiu)e'J<l>oTdr  (1.15) 

Di 


regardless  of  the  size  of  Tf  . 

We  know  that  given  a  certain  temporal  sampling  rate,  Tf  will  determine  the  frequency 
sampling  rate  or  discrete  Fourier  transform  (DFT)  binwidth,  i.e., 


(1.16) 


for  oj  in  radians/second.  So  we  can  see  that  the  acoustic  propagation  environment  (charac¬ 
terized  by  Ls  and  Ds )  will  dictate  the  maximum  coarseness  of  the  (w,/?)  search  space  and 
as  a  result  influence  our  calculation  of  the  steering  vectors  of  equations  1.3  and  1.4. 
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1.2.1. 5  Broadband  Definition 


Analogous  with  the  definitions  given  in  Wax,  Shan,  and  Kailath  (1984)  and  Hudson 
(1981),  we  define  a  source  to  be  narrowband  with  respect  to  some  array/source  configura¬ 
tion  if 


B« 
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Ls-Df 


where  Ls  and  Ds  are  from  section  1.2. 1.4.  It  is  broadband  otherwise. 


(1.17) 


1.2.2  Time-Variant  Channel  Filter  Model 

Realistically,  however,  the  channel  filters  rarely  are  time  invariant  (we  have  included 
the  time-variant  impulse  responses  in  the  (  )s  in  figure  1.1).  There  might  be  random 
perturbations  in  the  channel  caused  by  moving  scatterers,  moving  sea  surface,  or  a  mov¬ 
ing  source  (Knight,  Pridham,  and  Kay,  1981).  In  this  case,  x(t;fisi)  and  yt  ( t )  will  not  be 
WSS  (Weinstein,  1978)  since  x(t;fisi)  will  not  be  the  linear  convolution  of  j(t)  and  h(t;fisi ) 
and  so  the  concept  of  a  CSDM  then  has  no  meaning.  To  deal  with  this  situation,  we 
usually  assume  that  a  segment  of  data  is  a  windowed  portion  of  WSS  data  and  we  try  to 
estimate  the  statistics  of  this  “imaginary”  data  using  the  segment  as  in  equation  1.8.  The 
more  slowly  varying  the  filters,  the  longer  we  can  make  the  window  and  the  better  our 
estimates  get. 

Problems  can  result  when  K  in  equation  1.8  is  large  compared  to  the  time-varying 
nature  of  the  propagation  environment.  Usually  a  lack  of  coherence  between  sensors 
results,  which  manifests  itself  in  the  calculation  of  R{w>),  e.g.,  for  the  environment  of 
figure  1-1  we  will  not  get  the  cross  sensor  statistics  of  equation  1.2.  Usually  for  large 
arrays  and  long-source  ranges,  the  random  channel  perturbations  cause  this  lack  of  coher¬ 
ence  (Morgan  and  Smith,  1990).  For  shorter  ranges,  the  moving  source  becomes  the 
primary  cause  (Gerlach,  1978,  and  Patzewitsch  and  Srinath,  1978).  In  this  case,  we  would 
expect  a  certain  degree  of  performance  degradation  from  the  spatial  processors  of  equa¬ 
tions  1.3  and  1.4  resulting  in  a  loss  of  output  signal-to-noise  ratio  (SNR),  an  increase  in 
peak  location  variance,  or  bias  etc.  Indeed  if  the  time-varying  nature  of  the  channel  is 
very  severe,  we  might  not  be  able  to  do  any  spatial  spectral  estimation  via  the  methods  of 
equations  1.3,  1.4,  and  1.8. 

However,  attempts  have  been  made  at  using  this  time-varying  nature  as  a  target  dis¬ 
criminant,  e.g.,  the  moving  source  situation  in  Weinstein  (1978).  Viewing  the  problem  as 
in  figure  1-1  might  allow  us  also  to  use  time-variant-filter  theory  for  source  localization 
(Chen,  1984,  and  Kitagawa  and  Gersch,  1985). 


1.3  CORRELATION  BETWEEN  FOURIER  COEFFICIENTS 

In  our  calculations  of  the  CSDM,  we  only  have  determined  correlations  between 
Fourier  coefficients  of  the  same  frequency.  Since  the  signal  has  some  frequency  extent, 
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one  might  wonder  if  there  are  correlations  between  Fourier  coefficients  that  we  might  use 
to  enhance  array  processing  performance.  It  can  be  shown  (Papoulis,  1984;  Blachman, 
1957;  and  Hodgkiss  and  Nolte,  1976)  that  if  we  define  xt  ( (om )  and  as  in  equation 

1.6  (without  the  noise  component),  then 

lira  4  x'((on)}  =  0  V  i,k,  and  m  *  n  (1.18) 

T  —  =s  / 

where  m  and  n  are  frequency  indexes.  That  is,  the  “cross-frequency”  power  spectrum  of  a 
WSS  process  is  zero,  i.e.,  the  Fourier  coefficients  at  different  frequencies  are  uncorre¬ 
lated.  However,  these  coefficients  in  fact  are  correlated  for  finite  T.  An  expression  of  the 
correlation  easily  can  be  derived  using  equation  1.2  and  equation  1.6  (Blachman,  1957, 
and  Hodgkiss  and  Nolte,  1976)  giving 

T  T 

E{x.(com)xmk((on}  =  J  J  rXjk(t-s)  e10>mt  JMnSdtds 

T  T 

2  2  (1.19) 

■  h  / q W«(tt>>sinc[(ai  "2n,m)r}/'i{(ai 
-  00 

where  r  (f  -  s)  is  the  cross  correlation  function  of  the  ith  and  kth  processes.  We  can  see 
that  equation  1.19  approaches  equation  1.18  when  T  -»  °°  for  m  *  n  and  approaches 
equation  1.5  for  m  =  n. 

The  conclusions  reached  in  Hodgkiss  and  Nolte  (1976)  indicate  that  if  a2psix^ski(Ji)  is 

fairly  constant  over  the  interval  where  sinc[— -”'r‘)r]sinc[ has  appreciable  value, 
then  equation  1.19  approximates  the  correlation  value  for  T  -*  This  result  is  useful  for 
certain  frequency-domain  array  processing  implementations  (Hodgkiss,  1979).  Since  this 
assumption  obviously  is  dependent  on  the  propagation  environment,  it  is  unclear  how 
applicable  it  is  to  non-plane-wave  conditions.  A  certain  degree  of  correlation  actually 
might  enhance  array  processing  performance  if  used  correctly,  e.g.,  in  a  “cross¬ 
frequency”  power  spectral  density  matrix. 

1.3.1  Frequency  Averaging 

Interestingly,  frequency  averaging  as  in  equation  1.8  appears  to  be  a  faulty  method  of 
estimating  equation  1.19.  To  see  this,  we  define  the  continuous  frequency  average  as 

A 

*>*)>  =  Jim,  JX  /  *,70KKfoK)d*o  (1-2°) 

-A 

where 
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(121) 


*J“>J  -  I  *M'>e''“m'de 


*1  (°0  = 


J  xk(oi)  e’con0idoi 


(1.22) 


~7+'° 


Then 


/»  2  2 


£/AK)i;k)>=  Hm  T7  /  J  J^+jo)jc*(5+*o>c"71 


jMo>m-o>n)e-jo>n,'e  jo>ns ^  ^  ^  (123) 


where  t  =  €  -  tn  and  s  =  a  -  t  .  Interchanging  the  order  of  integration,  we  Ret 


T  jA~*x  _A 


Now  if,  a>  =  0)  and  using  the  ergodicity  property  of  x  ( t ),  i.e., 

m  n  1 » * 


lim  77 

,4  -*•  00  Z/i 


/H 

T7  /  -t  ('  + + I„)*3  =  -J)  • 


(1.25) 


we  get 


£/a^KKK)>=  J  I  rx.k(t-s)fM,'s)dtds  . 


(1.26) 


Then  (Thomas,  1984) 


y£{x  K)t;(«  )}  =  ^W^C")  •  (1-27) 


This  is  the  familiar  result  from  section  1.2. 1.3.  Now  if,  a)m  *  co  ,  we  want  to  arrive  at 
equation  1.19.  However,  in  general 

1  a 

A™*  2A  I  *,('  +  ^  +  ^'"'“"'""’<*0  *  r,„('-s)  •  (1.28) 

This  is  easy  to  see  since  rXjk(t  -  s )  is  real  and  equation  1.28  is  a  windowed  Fourier 
transform  of  x At  +  t  )x  (s  +  t  )  over  the  frequency  variable  (o  -c o  .  Since  x.(t  +  :  )x 
(s  +tQ)  is  not  necessarily  an  even  function,  the  left  side  of  equation  1.28  in  general  is 
complex. 

So  it  appears  that  frequency  averaging  cannot  be  used  to  estimate  correlation  between 
Fourier  coefficients  at  different  frequencies.  Equation  1.19  indicates  one  approach  is  to 
take  a  two-dimensional  Fourier  transform  of  an  estimate  of  the  correlation  function. 
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2  0  METHODS  OF  BROADBAND  ARRAY  PROCESSING 


In  this  section,  we  introduce  approaches  to  recently  developed  broadband  array  proc¬ 
essing  that  take  advantage  of  the  frequency  extent  of  the  source  process  to  estimate  its 
spatial  spectrum.  The  frequency-domain  methods  essentially  collapse  a  region  of  the  fre¬ 
quency  axis  from  the  (o),/I)  search  space  and  then  search  over  the  (fi)  region,  while  the 
time-domain  methods  use  the  linear  systems  model  of  the  propagation  environment  to 
determine  matched  filters  for  the  spatial  searc.t  space.  Mainly,  we  will  be  concerned  with 
MVDR  processing. 

2.1  BROADBAND  MODEL 

We  use  the  following  standard  model  of  the  cross  spectral  density  matrix  of  equation 

R(co)  =  Aipj)  Ps(co)  AB(co)  +  Q(co)  .  (2.1) 

If  there  are  N  sources,  then  A  (to)  is  an  M  x  N  matrix  whose  columns  are  the  signal 
vectors  of  the  sources.  Pj(co)  is  an  N  x  N  matrix  whose  diagonal  elements  are  the  power 
spectral  densities  of  each  source  and  off-diagonal  elements  are  the  respective  cross-source 
spectral  densities.  If  the  sources  are  mutually  uncorrelated,  Ps (w)  is  a  diagonal  matrix. 
(2(to)  is  the  noise  matrix. 

2.2  INCOHERENT  FREQUENCY-DOMAIN  ARRAY  PROCESSING 

The  simplest  and  most  common  approach  to  broadband  array  processing  is  to  decom¬ 
pose  the  process (es)  received  at  the  sensors  into  narrowband  components  via  the  Fourier 
transform,  estimate  the  spatial  spectrum  using  equations  1.3,  1.4,  and  1.5,  and  then  inco¬ 
herently  sum  up  the  spectrums  across  the  frequency  range  of  interest  (Tolstoy  and  Porter, 
1986),  i.e., 

Pi_  conv(#)  =  2PConM^  (2.2) 

/ 

or 

Pijnvdffl  =  X^mvdr  .  (2.3) 

i 

The  success  of  this  method  rests  on  the  fact  that  spurious  peaks  resulting  from  sidelobes 
incoherently  add  across  frequency,  while  actual  peaks  resulting  from  the  sources  coher¬ 
ently  add.  Thus,  we  might  be  able  to  effectively  process  broadband  signals  with  airays 
whose  “half-wavelength”  sensor  spacing  is  much  greater  than  the  wavelengths  associated 
with  the  spectral  composition  of  the  broadband  process.  This  allows  us  to  use  sparse 
arrays  and  still  avoid,  to  a  certain  extent,  spatial  aliasing. 
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2.2.1  Simulations 


Initially,  our  simulations  consist  of  two  mutually  uncorrelated  sources  impinging  upon 
a  16-element  line  array  with  7.5-meter  sensor  spacing,  i.e.,  the  half  wavelength  spacing 
for  100  Hz  assuming  1500  meter/second  sound  speed.  Both  sources  are  plane  waves  with 
0°  declination  angle.  Source  1  emanates  from  180°  azimuth  and  source  2  from  186°  azi¬ 
muth,  where  180°  represents  broadside.  Both  sources  are  bandpass,  Gaussian  processes 
with  bandwidth  B  =  40  Hz  with  100-Hz  center  frequency.  We  also  include  sensor  noise 
(temporally  and  spatially  uncorrelated).  The  power  spectral  density  values  of  the  source 
and  noise  processes  are  0  dB  within  the  80-  to  120-Hz  band.  We  note  that  the  sources  are 
plane  waves  for  convenience  only  and  the  results  given  here  are  applicable  to  more  gen¬ 
eral  propagation  environments. 

We  spectrally  decomposed  the  data  into  0.934-Hz  binwidths  and  constructed  the  nar¬ 
rowband  CSDMs  with  56  independent  fast  Fourier  transforms  (FFTs),  which  comprised 
approximately  1  minute  of  data.  Figures  2-1  and  2-2  are  the  results  of  the  narrowband 
MVDR  processor  for  100  Hz  and  120  Hz  respectively  for  10  runs,  where  each  run  repre¬ 
sents  an  estimate  of  the  spatial  spectrum  of  the  signal  and  noise  processes  for  some 
sample  function.  For  both  frequencies  we  see  the  variance  in  the  array  response  in  the 
direction  of  the  signals  and  most  notably  in  the  nonsignal  directions.  This  variance  is  due 
to  inaccuracies  in  estimating  the  CSDM. 

We  get  the  results  shown  in  figures  2-3  and  2-4  by  incoherently  summing  the  array 
responses  across  the  80-  to  120-Hz  frequency  range,  as  in  equations  2.2  and  2.3  respec¬ 
tively,  for  the  10  independent  runs.  We  first  note  that  the  conventional  processor  does  not 
quite  resolve  the  two  sources,  while  the  MVDR  processor  clearly  does.  We  also  see  a 
dramatic  decrease  in  MVDR  array  response  variance  over  the  results  shown  in  figures  2-1 
and  2-2.  Certainly  one  of  the  advantages  of  incoherent  processing  might  be  the  enhanced 
detection  of  low-SNR  broadband  signals  in  the  presence  of  strong  narrowband  or  broad¬ 
band  interferers  by  this  decrease  in  variance  of  the  side  regions  of  the  array  response. 

2.3  COHERENT  FREQUENCY-DOMAIN  ARRAY  PROCESSING 

Many  methods  have  been  proposed  recently  that  attempt  to  effectively  use  the  broad¬ 
band  nature  of  the  source  to  enhance  array  processing  performance  (Wang  and  Kaveh, 
1985;  Krolik  and  Swingler,  1989;  Buckley,  1987;  Buckley  and  Griffiths,  1988;  and  Nawab, 
Dowla,  and  Lacoss,  1985).  Here  we  will  concentrate  on  the  MVDR-based  broadband 
method  presented  in  Krolik  and  Swingler  (1989),  which  is  essentially  derived  from  the 
approach  given  in  Wang  and  Kaveh  (1985).  In  this  method,  the  narrowband  CSDMs  are 
estimated  as  in  the  incoherent  technique  of  section  2.2.  Then  they  are  presteered  to  some 
direction  of  look  and  coherently  summed  across  the  desired  spectrum.  The  resulting 
steered  matrix  is  then  used  in  an  MVDR  algorithm. 
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Figure  2-1.  Narrowband  MVDR  processor  results  at  100  Hz  with  56  independent  FFTs 
per  CSDM. 
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Figure  2-2.  Narrowband  MVDR  processor  results  at  120  Hz  with  56  independent  FFTs 
per  CSDM. 
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Figure  2-3.  Incoherent  conventional  processor  results  with  56  independent  FFTs 
per  CSDM. 
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Figure  2-4.  Incoherent  MVDR  processor  results  with  56  independent  FFTs 
per  CSDM. 
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To  begin,  we  expand  equation  2.2  using  equation  1.3,  i.e., 

pL conv^)  =  ^  JaeH{a)l-p)R(oj)e{a):,p)  . 

Now, 

*(<».; £)  =  E(a)-P)l  , 
where  E(aj;P)  is  the  diagonal  matrix 


e2(a>rP) 

0 


0 


(2.4) 


(2.5) 


(2.6) 


and  1  is  an  M  x  1  vector  of  ones.  So, 


or 


(2.7) 


(2.8) 


where  Rst(a>:,P)  is  the  “steered"  CSDM  (SCSDM).  We  can  see  that  if  the  signal  vector  of 
some  source  described  by  a  column  of  A(co.)  of  equation  2.1  is  modeled  by  e(w:,P),  then 
this  signal  is  mapped  onto  the  1  vector  in  the  operation  of  equation  2.7.  Then,  these 
“mapped”  components  of  the  broadband  signal  are  coherently  added  across  frequency  in 
equation  2.8  while  the  other  sources  represented  in  A(w.)  and  the  noise  components  of 
Q(<y.)  are  incoherently  summed. 


Now,  recognizing  the  similarity  between  equation  2.8  and  equation  1.3,  we  can  see 
that  an  MVDR-based  method  or  coherent  MVDR  technique  that  uses  the  SCSDMs  results 
in 


Pc  mvdr($)  “ 


(2.9) 


We  note  that  the  interpretation  of  equation  2.2  given  by  equation  2.8  implies  that 
Pi_com(P)  =  Pc_conv(&  while  lt  is  obvious  that  Pc  mvdr(fi)  is  markedly  different  from 
Pj  mvdr(/f)  •  We  define  ^  Est{o).\  fi)  as  the  coherent  CSDM  or  Rc(fi) .  Then 
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1 


(2.10) 


c  mvdr 


<0  = 


i  TR-:m 


An  obvious  drawback  of  the  method  described  by  equation  2.9  is  that  SCSDMs  must 
be  estimated,  summed,  and  inverted  for  each  direction  of  look.  This  might  be  a  consider¬ 
able  computational  burden  for  large  arrays.  However,  the  main  benefit  of  this  method 
over  the  incoherent  technique  of  equation  2.3  is  that  there  are  more  statistical  degrees  of 
freedom  available  to  estimate  Rc(fi)  than  are  required  to  estimate  each  narrowband 
CSDM  (Krolik  and  Swingler,  1989).  So,  given  some  Fourier  transform  length,  we  might 
be  able  to  estimate  a  “statistically  sufficient”  R  (fi)  using  fewer  Fourier  transform  updates 
than  the  number  required  to  estimate  the  narrowband  CSDMs.  This  would  be  of  value  in 
highly  nonstationary  environments  as  described  in  section  1.2.2.  Also,  coherent  MVDR 
easily  can  be  implemented  in  a  parallel  processing  environment  since  the  coherent 
CSDMs  for  all  directions  of  look  can  be  constructed  and  inverted  independently  of  one 
another. 


We  also  note  that  an  appropriate  steering  matrix  for  broadband  coherent  MVDR 
matched-field  array  processing  might  be  one  whose  elements  are  inverses  of  the  respective 
channel  frequency  response  values  rather  than  simply  conjugates  as  described  in  section 
1.0.  For  plane-wave  processing,  these  two  cases  are  identical. 


2.3.1  Simulations 

We  now  analyze  the  synthetic  broadband  data  described  in  section  2.2.1  using  the 
coherent  broadband  MVDR  technique  of  equation  2.9.  Spectrally  decomposing  the  data  as 
before,  we  get  the  results  of  figure  2-5  using  56  independent  FFTs  per  CSDM.  We  note 
the  similarity  between  figure  2-5  and  the  incoherent  MVDR  results  of  figure  2-4.  We  can 
conclude  that  nothing  is  gained  from  the  increased  computational  expense  of  coherent 
MVDR  if  enough  data  are  available  to  compute  well-conditioned  narrowband  CSDMs 
(Giannella  and  Schultheiss,  1990). 

Next,  we  construct  the  matrices  using  only  one  outer  product  or  approximately  1  sec¬ 
ond  of  data.  Technically,  these  matrices  are  not  CSDM  estimates  since  no  frequency 
averaging  is  done.  Since  the  number  of  frequency  bins  is  much  greater  than  16,  Rc( ft)  is 
full  rank  and  we  get  the  results  shown  in  figure  2-6.  We  note  the  considerable  increase  in 
variance  of  the  power  in  the  side  regions  as  well  as  in  the  direction  of  the  two  signals.  We 
also  see  an  increase  in  peak  bias  in  these  two  directions.  However,  both  signals  are  still 
well  localized.  Figures  2-7  and  2-8  represent  results  for  two  and  three  FFTs  per  “CSDM” 
respectively.  We  see  the  gradual  decrease  in  variance  and  peak  bias.  In  fact,  the  results  of 
figure  2-8  appear  similar  in  quality  to  the  100-Hz  narrowband  result  of  figure  2-1.  There 
are  no  incoherent  MVDR  results  that  we  can  compare  to  figures  2-6,  2-7,  and  2-8  since 
the  narrowband  CSDMs  are  singular  in  these  cases. 
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Figure  2-5.  Coherent  MVDR  processor  results  with  56  independent  FFTs 
per  CSDM. 
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Figure  2-6.  Coherent  MVDR  processor  results  with  one  FFT  per  “CSDM.” 
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Figure  2-7.  Coherent  MVDR  processor  results  with  two  independent  FFTs 
per  “CSDM.” 
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Figure  2-8.  Coherent  MVDR  processor  results  with  three  independent  FFTs 
per  “CSDM.” 
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2.3.2  Correlated  Sources 


An  additional  interesting  aspect  of  the  coherent  MVDR  algorithm  is  that  it  is  able  to 
localize  completely  correlated  broadband  sources  as  long  ns  the  signal  vectors  are  not 
equal  or  simple  scaler  multiples  of  one  another.  This  is  essentially  due  to  the  fact  that 
presteering  and  frequency  averaging  results  in  coherent  addition  only  for  sources  modeled 
by  the  steering  matrix  regardless  of  the  power  spectral  density  of  the  source.  The 

greater  the  spectral  extent  of  the  averaging,  the  better  the  performance.  See  Wang  and 
Kaveh  (1985)  for  a  more  general  explanation. 

2.3.2. 1  Simulations 

For  simulation  purposes  we  modify  the  source  characteristics  of  section  2.2.1  so  that 
the  two  sources  are  completely  correlated.  Figure  2-9  is  the  result  of  incoherent  MVDR 
for  this  case.  Interestingly,  two  sources  are  apparently  detected.  However  this  is  due  to  the 
fact  that,  because  of  the  regularity  of  the  array,  the  resulting  correlated  signal  vector  (the 
sum  of  the  two  individual  signal  vectors)  has  two  prominent  sidelobes  in  these  directions. 
Note  that  there  is  considerable  peak  bias  for  all  the  trials.  Conversely,  figure  2-10  results 
from  coherent  MVDR.  While  there  is  a  considerable  drop  in  output  power,  there  is  no 
significant  peak  bias. 

2.3.3  Multiple  Presteering 

The  coherent  MVDR  method  requires  a  matrix  inverse  for  each  direction  of  look  since 
the  narrowband  CSDMs  are  presteered  to  only  one  direction.  A  modification  of  this  tech¬ 
nique  would  presteer  for  multiple  look  directions. 

Let  B(q).)  be  an  M  x  P  matrix  where  each  of  the  P  columns  is  a  steering  vector 
representing  some  direction  of  look.  We  want  to  find  a  steering  matrix  T((a>.;  l?(a>.))  such 
that 


Tsi((o-B(cd))B((d)  =  G,  (2.11) 

where  G  is  an  M  x  P  matrix  whose  columns  represent  vectors  upon  which  we  want  the 
respective  columns  of  B{a).)  to  be  mapped.  Relating  equation  2.5  to  equation  2.11  for  the 
P  =  1  case,  we  see  that  Tst(w:,B( a>.))-  E  '(ca.;/?),  fi(o>. )  =  e(a).\P),  and  G=  1.  However,  P 
can  be  greater  than  1  and  it  is  beneficial  to  have  P  as  large  as  possible  to  decrease  the 
number  of  matrix  inversions.  In  fact,  if  we  choose  P  =  M  we  get 

TJu-Bio)'))  =  GB\w .)  (2.12) 

if  B{a)t)  is  full  rank.  Note  that  when  implementing  equation  2.12,  we  might  encounter 
numerical  difficulties  if  our  presteering  directions  are  spatially  very  close  since  the 
inversion  of  B(wt)  might  become  unstable,  i.e.,  fl(cy.)  might  be  numerically  singular. 
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Figure  2-9.  Incoherent  MVDR  processor  results  for  two  perfectly  correlated 
sources. 
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Figure  2-10.  Coherent  MVDR  processor  results  for  two  perfectly  correlated  sources. 
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If  1  <  P  <  M  we  could  use  the  pseudo-inverse  of  B( a»; )  to  calculate  fi(ru.)) t  If  P  >  M 

there  is  no  solution  to  this  overdetermined  system  of  equations.  However,  we  might  still 
obtain  adequate  presteering  if  we  use  the  least  squares  solution  to  the  problem  as 

r„(©j; «(<».))  =  GBV/)|fl(<u,)«V,)r'  ■  (2-13) 


How  well  we  can  presteer  and  coherently  add  across  frequency  if  P  »  M  is  unclear. 

Now,  using  equation  2.12  or  equation  2.13  for  our  steering  matrix,  we  get  the 
multiply-steered,  coherent  CSDM 


R(B) 


X  7 «(©.;  B(co))R(oj)Tt((o.-  B(co .))  . 


(2.14) 


Then, 


c  mvdrV£>fc 


(«,)  - 


(2.15) 


where  g  is  the  kth  column  of  G.  So  for  each  matrix  inverse  we  can  steer  over  P  directions 
where  1  <  k  <  P.  The  multiply-steered  coherent  MVDR  method  is  then  P  times  less 
computationally  intensive  than  the  technique  of  equation  2.9.  So,  generally,  the  larger  the 
number  of  sensors,  the  more  directions  we  can  presteer  towards  in  one  RJJP) .  The  compu¬ 
tational  overhead  of  coherent  MVDR  might  not  be  severe  for  large  arrays. 

We  have  complete  freedom  in  choosing  G.  A  convenient  form  would  be  an  orthogonal 
or  orthonormal  matrix  if  P  <  M.  Then  if  Tsl(o).\  Z?(<w.))  contains  more  than  one  source  signal 
vector  from  A(a).)  of  equation  2.1,  Rc(Jf)  would  contain  orthogonal  components.  This  in 
turn  could  enhance  nulling  performance.  We  could  also  use  a  G  matrix  with  M  orthogonal 
basis  columns  and  P  -  M  columns  spanned  by  the  basis  vectors  if  P  >  M. 

We  could  also  use  multiple  presteering  for  coherent  conventional  processing.  Then 
multiply-steered  coherent  conventional  processing  would  not  equal  incoherent  conven¬ 
tional  results.  It  is  unclear  what  the  performance  differences  would  be  between  single  and 
multiple  direction  coherent  MVDR  and  conventional  processing.  Certainly,  presteering 
and  coherently  summing  two  or  more  signal  vectors  in  one  coherent  CSDM  will  affect 
broadband  array  processing  performance. 


2.3.4  Coherent  Processing  Applied  to  MUSIC 

We  have  confined  our  attention  to  coherent  MVDR  mainly  for  instructional  purposes. 
However,  it  is  a  simple  matter  to  use  the  coherent  CSDM  Rc(fi)  from  section  2.3  in  the 
MUSIC  algorithm.  If  we  compute  the  eigen  decomposition  of  Rc(fi)  as 

M 

m  -  z  k*<&  •  <216) 

i  =  I 
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where  vci  are  the  eigenvectors  and  Xa the  associated  eigenvalues  with  Acl  >  Xc2  >  ...  >  XcM, 
then  the  well-known  MUSIC  algorithm  produces 


Pc  music  (P)  = 


M  -  n 


y  irv  v"i 

ci  ci 


i  =  1 


(2.17) 


where  n  is  the  user-specified  order.  If  we  were  able  to  model  the  noise  matrix  component 
of  R(fi) ,  we  could  also  implement  a  generalized  MUSIC  algorithm  as  in  Wang  and  Kaveh 
(1985).  Foi  multiply-steered  CSDMs  we  get 


P  (e  ) 

c  music'oA:'' 


M-n 


1  A 


i  =  i 


(2.18) 


unfortunately,  order  selection  in  MUSIC  can  be  quite  problematic. 
2.3.4. 1  Simulations 


We  examine  the  narrowband  “CSDMs”  constructed  from  three  independent  FFTs  per 
matrix  used  to  produce  the  results  of  figure  2-8.  Figures  2-11,  2-12,  and  2-13  are  the 
MUSIC  results  from  equation  2.17  for  orders  one,  two,  and  three  respectively.  Since  we 
are  incoherently  summing  signal  components  across  frequency  in  RJJf)  when  we  are  not 
steered  towards  a  component,  one  might  assume  an  order  of  one  would  be  appropriate 
because  there  should  be  at  most  one  dominant  correlated  component  in  Rc(ft) .  Figure  2-1 1 
indicates  that  this  is  not  the  case  and  an  order  of  two  is  still  required  since  spectral 
averaging  is  not  sufficient  to  completely  “whiten”  a  cor  "elated  component.  Comparison  of 
figure  2-12  with  figure  2-8  reveals  that  coherent  MUSIC  offers  greater  resolution  in  the 
signal  directions  and,  interestingly,  less  variance  in  the  side  regions. 


2.4  TIME-DOMAIN  ARRAY  PROCESSING 

Methods  have  been  proposed  recently  that  localize  acoustic  sources  using  time-domain 
techniques  (Clay,  1987;  Li  and  Clay,  1987;  Clay  and  Li,  1988;  and  Hodgkiss  and  Brienzo, 
1990).  These  methods  use  matched-filter  or  inverse-filter  approaches  to  essentially 
remove,  to  a  certain  extent,  the  propagation  environment  from  the  signals  received  at 
each  of  the  sensors  and  then  to  use  various  processing  methods  such  as  estimating  corre¬ 
lation  statistics  to  localize  sources. 

Figure  2-14  is  a  generalized  diagram  of  these  techniques.  y .  (t)  represents  the  data 
received  at  the  ilh  sensor  (see  figure  1-1),  and  h(t;(fsi)  is  the  matched  or  inverse  filter 
impulse  response  for  the  ilh  sensor  focused  or  steered  to  a  spatial  location  specified  by  fisi. 
We  have  also  included  a  time-varying  steering  filter  h,(t;fisil)  in  the  diagram  to  indicate 
the  possible  extension  to  a  spatial  and/or  temporal  nonstationary  environment. 
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Figure  2-11.  Coherent  MUSIC  results  with  n  =  1  and  three  independent  FFTs 
per  “CSDM.” 
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Figure  2-12.  Coherent  MUSIC  results  with  n  =  2  and  three  independent  FFTs 
per  “CSDM.” 


22 


OUTPUT  POWER  (dB) 


135  140  145  150  155  160  165  170  17c  1E0  185  190  195  200  205  210  215  220  225 

AZIMUTH 

Figure  2-13.  Coherent  MUSIC  results  with  n  =  3  and  three  independent  FFTs 
per  “CSDM.” 
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Figure  2-14.  Generalized  time-domain  array  processor  diagram. 
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To  compute  the  focusing  filter,  the  channel  impulse  response  is  first  required.  It  can 
be  estimated  analytically  with  time-domain  propagation  models  or  via  Fourier  synthesis 
using  frequency-domain  propagation  techniques  (Hodgkiss  and  Brienzo,  1990),  or,  if  pos¬ 
sible,  determined  experimentally  using  an  impulsive  source  (Clay,  1987;  Li  and  Clay, 
1987;  Clay  and  Li,  1988)  The  outputs  of  these  filters  are  then  processed  using  simple 
summing  for  conventional  processing  as  in  Hodgkiss  and  Brienzo  (1990),  or  cross  corre¬ 
lating  between  sensors  as  in  Clay  (1987),  Li  and  Clay  (1987),  and  Clay  and  Li  (1988). 
Future  work  could  include  a  tapped-delay-line-based  MVDR  technique  as  in  Frost  (1972), 
which  is  also  useful  in  highly  nonstationary  environments.  In  fact  there  are  similarities 
between  a  Frost-based,  time-domain  matched-field  MVDR  processor  and  the  coherent 
MVDR  technique.  The  output  of  the  processor  can  then  be  sent  to  a  postprocessor  that 
might  spectrally  decompose  the  data,  compute  energy  (Hodgkiss  and  Brienzo,  1990),  or 
do  spectral  peak  detection  (Clay,  1987;  Li  and  Clay,  1987;  Clay  and  Li,  1988). 

2.4.1  Matched-Filter  Method 

In  the  matched-filter  approach,  each  of  the  steering  filters  of  figure  2-11  has  an 
impulse  response  that  is  the  time  reverse  of  the  respective  channel  impulse  response  with 
an  associated  time  shift  to  ensure  causality.  The  data  received  at  each  sensor  are  con¬ 
volved  with  the  impulse  response  of  the  matched  filter.  While  this  method  will  not  entirely 
remove  the  environment  from  the  signal,  the  more  complicated  the  environment  (the 
longer  the  channel  impulse  response),  the  better  the  performance  of  this  approach  (Clay, 
1987).  It  is  also  simple  and  very  robust. 

2.4.2  Inverse-Filter  Method 

A  time-domain  analogue  to  the  frequency-domain  techniques  of  equations  1.3  and  1.4 
is  the  inverse-filter  method.  Here  the  steering  filters  are  computed  such  that 

=  (5(0,  (2.19) 

where  *  denotes  linear  convolution  and  8(0  is  the  Dirac  delta  function.  So  if  we  steer 
directly  at  a  signal  s(t)  with  power  spectral  density  o2s{(S),  the  power  spectral  density  due 
to  the  signal  at  each  focusing  filter  output  also  will  be  a?(w)- 

As  noted  in  Clay  and  Li  (1988)  for  the  discrete-time  case,  care  must  be  taken  if  the 
channel  transfer  function  is  not  minimum  phase,  i.e.,  the  zeros  of  the  Laplace  transform 
of  h(t\Psi )  are  not  all  in  the  left  half  of  the  complex  plane,  or,  for  the  discrete-time  case, 
the  zeros  of  the  z-transform  of  h(n;Psi)are  not  all  inside  the  unit  circle.  The  causal  inverse 
filter  would  then  be  unstable,  i.e.,  the  impulse  response  would  not  be  absolutely 
integrable  (summable).  Standard  zero  reflection  .and  appropriate  time-shifting  can  solve 
this  problem,  e.g.,  see  Oppenheim  and  Schafer  (1989)  for  the  discrete-time  case.  This 
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technique  obviously  requires  more  processing  than  the  matched-filter  method.  However,  if 
the  channel  transfer  function  has  a  zero  on  the  imaginary  axis  (the  unit  circle  for  the 
discrete-time  case),  no  stable  inverse  filter  can  be  found.  An  interesting  topic  would  be  to 
determine  analytically  the  conditions,  e.g.,  the  boundary  conditions  on  the  wave  equation, 
that  result  in  a  minimum-phase  channel. 

We  can  now  recognize  a  benefit  of  the  frequency-domain  coherent  MVDR  algorithm 
over  a  time -domain,  Frost-based  MVDR  method.  Since  we  are  interested  only  in  calculat¬ 
ing  a  particular  frequency  response  in  equation  2.7  over  the  coi  of  interest  rather  than  in 
calculating  an  impulse  response,  we  do  not  have  to  worry  about  causality  or  stability  of 
the  inverse  filter. 
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